Plot models
Max Lindmark
04 January, 2023
# Load libraries, install if needed
library(tidyverse)
#> Warning: package 'tidyr' was built under R version 4.0.5
library(tidylog)
library(RCurl)
library(RColorBrewer)
#> Warning: package 'RColorBrewer' was built under R version 4.0.5
library(sdmTMB)
library(devtools)
library(data.table)
library(here)
library(stringr)
library(boot)
library(sdmTMB)
library(patchwork)
# Source code for map plots
# change to url once we have the final one
source("/Users/maxlindmark/Dropbox/Max work/R/marine-litter/R/functions/map-plot.R")
#> Warning: package 'sf' was built under R version 4.0.5
# Source utm function
source_url("https://raw.githubusercontent.com/maxlindmark/marine-litter/main/R/functions/lon-lat-utm.R")
# Continuous colors
options(ggplot2.continuous.colour = "viridis")
# Discrete colors
scale_colour_discrete <- function(...) {
scale_colour_brewer(palette = "Dark2")
}
scale_fill_discrete <- function(...) {
scale_fill_brewer(palette = "Dark2")
}Read data
# Read and make data long so that I can for loop through all categories
biom <- read_csv("data/west_coast_litter_biomass.csv") %>%
mutate(year_f = as.factor(year),
quarter_f = as.factor(quarter),
depth_sc = (depth - mean(depth)) / sd(depth)) %>%
rename(plastic = plast) %>%
pivot_longer(c(plastic, sup, fishery),
names_to = "model", # change this to model from litter category for plotting
values_to = "density") %>%
filter(model %in% c("fishery", "plastic", "sup")) %>%
mutate(model = str_to_title(model))
#> Rows: 609 Columns: 23
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (4): haul_id, ansse_area, trend_area, var
#> dbl (18): diverse, naturmaterial, plast, metall, glas, gummi, quarter, st_no...
#> lgl (1): balse_area
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: new variable 'year_f' (factor) with 10 unique values and 0% NA
#>
#> new variable 'quarter_f' (factor) with 2 unique values and 0% NA
#>
#> new variable 'depth_sc' (double) with 125 unique values and 0% NA
#>
#> rename: renamed one variable (plastic)
#>
#> pivot_longer: reorganized (plastic, sup, fishery) into (model, density) [was 609x26, now 1827x25]
#>
#> filter: no rows removed
#>
#> mutate: changed 1,827 values (100%) of 'model' (0 new NA)
num <- read_csv("data/west_coast_litter_numbers.csv") %>%
mutate(year_f = as.factor(year),
quarter_f = as.factor(quarter),
depth_sc = (depth - mean(depth)) /sd(depth)) %>%
rename(plastic = plast) %>%
pivot_longer(c(plastic, sup, fishery),
names_to = "model", # change this to model from litter category for plotting
values_to = "density") %>%
mutate(number = density * swept_area_km2) %>%
filter(model %in% c("fishery", "plastic", "sup")) %>%
mutate(model = str_to_title(model))
#> Rows: 609 Columns: 23
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (4): haul_id, ansse_area, trend_area, var
#> dbl (18): diverse, naturmaterial, plast, metall, glas, gummi, quarter, st_no...
#> lgl (1): balse_area
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: new variable 'year_f' (factor) with 10 unique values and 0% NA
#>
#> new variable 'quarter_f' (factor) with 2 unique values and 0% NA
#>
#> new variable 'depth_sc' (double) with 125 unique values and 0% NA
#>
#> rename: renamed one variable (plastic)
#>
#> pivot_longer: reorganized (plastic, sup, fishery) into (model, density) [was 609x26, now 1827x25]
#>
#> mutate: new variable 'number' (double) with 125 unique values and 0% NA
#>
#> filter: no rows removed
#>
#> mutate: changed 1,827 values (100%) of 'model' (0 new NA)
binom <- read_csv("data/west_coast_litter_biomass.csv") %>%
mutate(year_f = as.factor(year),
quarter_f = as.factor(quarter),
depth_sc = (depth - mean(depth)) / sd(depth)) %>%
rename(other = diverse,
natural = naturmaterial,
glass = glas,
metal = metall,
rubber = gummi) %>%
pivot_longer(c(other, natural, glass, metal, rubber),
names_to = "model", # change this to model from litter category for plotting
values_to = "density") %>%
filter(!model %in% c("fishery", "plastic", "sup")) %>%
mutate(present = ifelse(density == 0, 0, 1)) %>%
mutate(model = str_to_title(model))
#> Rows: 609 Columns: 23
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (4): haul_id, ansse_area, trend_area, var
#> dbl (18): diverse, naturmaterial, plast, metall, glas, gummi, quarter, st_no...
#> lgl (1): balse_area
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: new variable 'year_f' (factor) with 10 unique values and 0% NA
#>
#> new variable 'quarter_f' (factor) with 2 unique values and 0% NA
#>
#> new variable 'depth_sc' (double) with 125 unique values and 0% NA
#>
#> rename: renamed 5 variables (other, natural, metal, glass, rubber)
#>
#> pivot_longer: reorganized (other, natural, metal, glass, rubber) into (model, density) [was 609x26, now 3045x23]
#>
#> filter: no rows removed
#>
#> mutate: new variable 'present' (double) with 2 unique values and 0% NA
#>
#> mutate: changed 3,045 values (100%) of 'model' (0 new NA)
# Load pred grid
# pred_grid_west <- read_csv("data/pred_grid_west.csv") %>%
# mutate(year_f = as.factor(year),
# quarter_f = as.factor(1),
# depth_sc = (depth - mean(num$depth)) / sd(num$depth)) %>%
# mutate(X = X*1000,
# Y = Y*1000)Read model files
# Example
# Read model files
file_paths <- list.files(path = here::here("output/models"),
pattern = ".rds",
full.names = TRUE)
file_names <- gsub(pattern = ".rds$",
replacement = "",
x = basename(file_paths))
for(i in 1:length(file_names)){
assign(file_names[i], readRDS(file_paths[i]))
}Read csv files with predictions, sims and coefficients
pred_binom <- read_csv("output/dat_pred_binom.csv") %>% mutate(model = str_to_title(model))
#> Rows: 212900 Columns: 11
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (1): model
#> dbl (10): X, Y, year, depth, year_f, quarter_f, depth_sc, est, est_non_rf, e...
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: changed 212,900 values (100%) of 'model' (0 new NA)
sim_binom <- read_csv("output/dat_sim_binom.csv") %>% mutate(model = str_to_title(model))
#> Rows: 50 Columns: 7
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (1): model
#> dbl (6): year, est, lwr, upr, log_est, se
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: changed 50 values (100%) of 'model' (0 new NA)
sims_binom <- read_csv("output/dat_sims_binom.csv") %>% mutate(model = str_to_title(model))
#> Rows: 25000 Columns: 4
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (1): model
#> dbl (3): year, .value, .iteration
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: changed 25,000 values (100%) of 'model' (0 new NA)
coef_binom <- read_csv("output/dat_coef_binom.csv") %>% mutate(model = str_to_title(model))
#> Rows: 10 Columns: 6
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (2): term, model
#> dbl (4): estimate, std.error, conf.low, conf.high
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: changed 10 values (100%) of 'model' (0 new NA)
pred_num <- read_csv("output/dat_pred_num.csv") %>% mutate(model = str_to_title(model))
#> Rows: 127740 Columns: 12
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (1): model
#> dbl (11): X, Y, year, depth, year_f, quarter_f, depth_sc, est, est_non_rf, e...
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: changed 127,740 values (100%) of 'model' (0 new NA)
sim_num <- read_csv("output/dat_sim_num.csv") %>% mutate(model = str_to_title(model))
#> Rows: 30 Columns: 7
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (1): model
#> dbl (6): year, est, lwr, upr, log_est, se
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: changed 30 values (100%) of 'model' (0 new NA)
sims_num <- read_csv("output/dat_sims_num.csv") %>% mutate(model = str_to_title(model))
#> Rows: 15000 Columns: 4
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (1): model
#> dbl (3): year, .value, .iteration
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: changed 15,000 values (100%) of 'model' (0 new NA)
coef_num <- read_csv("output/dat_coef_num.csv") %>% mutate(model = str_to_title(model))
#> Rows: 33 Columns: 6
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (2): term, model
#> dbl (4): estimate, std.error, conf.low, conf.high
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: changed 33 values (100%) of 'model' (0 new NA)
pred_biomass <- read_csv("output/dat_pred_biomass.csv") %>% mutate(model = str_to_title(model))
#> Rows: 127740 Columns: 12
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (1): model
#> dbl (11): X, Y, year, depth, year_f, quarter_f, depth_sc, est, est_non_rf, e...
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: changed 127,740 values (100%) of 'model' (0 new NA)
sim_biomass <- read_csv("output/dat_sim_biomass.csv") %>% mutate(model = str_to_title(model))
#> Rows: 30 Columns: 7
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (1): model
#> dbl (6): year, est, lwr, upr, log_est, se
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: changed 30 values (100%) of 'model' (0 new NA)
sims_biomass <- read_csv("output/dat_sims_biomass.csv") %>% mutate(model = str_to_title(model))
#> Rows: 15000 Columns: 4
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (1): model
#> dbl (3): year, .value, .iteration
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: changed 15,000 values (100%) of 'model' (0 new NA)
coef_biom <- read_csv("output/dat_coef_biomass.csv") %>% mutate(model = str_to_title(model))
#> Rows: 33 Columns: 6
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (2): term, model
#> dbl (4): estimate, std.error, conf.low, conf.high
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: changed 33 values (100%) of 'model' (0 new NA)pal <- brewer.pal(n = 8, name = "Paired")Plot depth coefficients
coefs <- bind_rows(coef_binom, coef_num, coef_biom) %>%
filter(term == "depth_sc") %>%
separate(model, sep = "_", into = c("Model", "Category")) %>%
mutate(Category = str_to_title(Category))
#> filter: removed 65 rows (86%), 11 rows remaining
#> mutate: changed 11 values (100%) of 'Category' (0 new NA)
ggplot(coefs, aes(y = estimate, x = Category, ymin = conf.low, ymax = conf.high, color = Model)) +
geom_point(size = 2, position = position_dodge(width = 0.2)) +
geom_errorbar(width = 0, position = position_dodge(width = 0.2)) +
geom_hline(yintercept = 0, size = 0.75, color = "gray", linetype = 2) +
scale_color_manual(values = pal[c(3, 2, 1)]) +
labs(y = "Estimate") +
theme_plot() +
coord_flip() +
NULL
ggsave("figures/coefs.pdf", units = "cm", width = 10, height = 20)Plot binomial predictions
dd <- binom %>%
group_by(model, year) %>%
summarise(mean_present = mean(present)) %>%
filter(!model %in% c("Glass", "Rubber"))
#> group_by: 2 grouping variables (model, year)
#> summarise: now 50 rows and 3 columns, one group variable remaining (model)
#> filter (grouped): removed 20 rows (40%), 30 rows remaining
p1 <- sims_binom %>%
filter(!model %in% c("Glass", "Rubber")) %>%
filter(.iteration < 26) %>%
ggplot() +
facet_wrap(~model, ncol = 3, scales = "free") +
geom_ribbon(data = sim_binom %>% filter(!model %in% c("Glass", "Rubber")),
aes(year, est, ymin = lwr, ymax = upr), alpha = 0.2) +
geom_line(aes(year, .value, group = .iteration), alpha = 0.2) +
geom_line(data = dd, aes(year, mean_present), color = pal[8], size = 0.75, linetype = 2) + # add data
geom_line(data = sim_binom %>% filter(!model %in% c("Glass", "Rubber")),
aes(year, est), size = 0.75) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 5)) +
theme_plot() +
theme(plot.margin = unit(c(0, 0.4, 0, 0.4), "cm")) +
xlab('Year') +
ylab('Presence prob. in tow') +
NULL
#> filter: removed 10,000 rows (40%), 15,000 rows remaining
#> filter: removed 14,250 rows (95%), 750 rows remaining
#> filter: removed 20 rows (40%), 30 rows remaining
#> filter: removed 20 rows (40%), 30 rows remaining
year_diff <- sims_binom %>%
filter(!model %in% c("Glass", "Rubber")) %>%
filter(year %in% c(min(year), max(year))) %>%
pivot_wider(names_from = year, values_from = .value) %>%
mutate("2022-2013" = `2022` - `2013`) %>%
pivot_longer(c(`2013`, `2022`, `2022-2013`)) %>%
ungroup()
#> filter: removed 10,000 rows (40%), 15,000 rows remaining
#> filter: removed 12,000 rows (80%), 3,000 rows remaining
#> pivot_wider: reorganized (year, .value) into (2013, 2022) [was 3000x4, now 1500x4]
#> mutate: new variable '2022-2013' (double) with 1,500 unique values and 0% NA
#> pivot_longer: reorganized (2013, 2022, 2022-2013) into (name, value) [was 1500x5, now 4500x4]
#> ungroup: no grouping variables
diff_sum <- year_diff %>%
filter(name == "2022-2013") %>%
group_by(model) %>%
summarise(diff_lwr = round(quantile(value, probs = 0.025), digits = 3),
diff_med = round(quantile(value, probs = 0.5), digits = 3),
diff_upr = round(quantile(value, probs = 0.975), digits = 3)) %>%
mutate(label = paste("2022-2013 = ", diff_med, " [", diff_lwr, ",", diff_upr, "]", sep = ""))
#> filter: removed 3,000 rows (67%), 1,500 rows remaining
#> group_by: one grouping variable (model)
#> summarise: now 3 rows and 4 columns, ungrouped
#> mutate: new variable 'label' (character) with 3 unique values and 0% NA
# Left join so that we can color "sig"
diff_sum <- diff_sum %>% mutate(sig = ifelse(diff_lwr > 0 & diff_upr > 0, "sig.", "not sig."),
sig = ifelse(diff_lwr < 0 & diff_upr < 0, "sig.", sig))
#> mutate: new variable 'sig' (character) with one unique value and 0% NA
year_diff <- left_join(year_diff, diff_sum)
#> Joining, by = "model"
#> left_join: added 5 columns (diff_lwr, diff_med, diff_upr, label, sig)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 4,500
#> > =======
#> > rows total 4,500
p2 <- ggplot(data = filter(year_diff, !name == "2022-2013"), aes(name, value)) +
geom_violin(alpha = 0.5, fill = "grey30", color = NA) +
geom_jitter(height = 0, width = 0.05, alpha = 0.2, color = "grey30", size = 0.1) +
geom_boxplot(fill = NA, width = 0.1, size = 0.5, color = "grey5",
outlier.shape = NA) +
theme_plot() +
theme(plot.margin = unit(c(0, 0.4, 0, 0.4), "cm")) +
facet_wrap(~ model, scales = "free", ncol = 3) +
xlab('Year') +
ylab('Presence prob. in tow') +
NULL
#> filter: removed 1,500 rows (33%), 3,000 rows remaining
p3 <- ggplot(data = filter(year_diff, name == "2022-2013"), aes(name, value, fill = sig)) +
geom_violin(alpha = 0.4, color = NA) +
geom_jitter(height = 0, width = 0.05, alpha = 0.3, color = "grey30", size = 0.2) +
geom_boxplot(fill = NA, width = 0.1, size = 0.5, color = "grey5",
outlier.shape = NA) +
scale_fill_manual(values = pal[c(2, 6)], name = "") +
theme_plot() +
theme(plot.margin = unit(c(0, 0.4, 0, 0.4), "cm")) +
facet_wrap(~ model, scales = "free", ncol = 3) +
geom_text(data = diff_sum, aes(label = label), inherit.aes = FALSE,
x = -Inf, y = Inf, size = 2.5, hjust = -0.05, vjust = 2) +
xlab('') +
ylab('Diff. in presence prob. per tow') +
geom_hline(yintercept = 0, linetype = 2, size = 0.4) +
NULL
#> filter: removed 3,000 rows (67%), 1,500 rows remaining
p1 / p2 / p3 + plot_annotation(tag_levels = "A")
ggsave("figures/binom_index.pdf", units = "cm", width = 20, height = 20)Plot Poisson predictions
# Space
pred_num
#> # A tibble: 127,740 × 12
#> X Y year depth year_f quart…¹ depth…² est est_n…³ est_rf omega_s
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 3.15e5 6.25e6 2013 228 2013 1 0.693 -1.38 -1.38 0 0
#> 2 3.17e5 6.25e6 2013 228 2013 1 0.693 -1.38 -1.38 0 0
#> 3 3.19e5 6.25e6 2013 227 2013 1 0.675 -1.38 -1.38 0 0
#> 4 3.21e5 6.25e6 2013 228. 2013 1 0.702 -1.38 -1.38 0 0
#> 5 3.23e5 6.25e6 2013 226 2013 1 0.656 -1.38 -1.38 0 0
#> 6 3.25e5 6.25e6 2013 229. 2013 1 0.717 -1.38 -1.38 0 0
#> 7 3.27e5 6.25e6 2013 229. 2013 1 0.711 -1.38 -1.38 0 0
#> 8 3.29e5 6.25e6 2013 228. 2013 1 0.686 -1.38 -1.38 0 0
#> 9 3.31e5 6.25e6 2013 229 2013 1 0.711 -1.38 -1.38 0 0
#> 10 3.33e5 6.25e6 2013 229 2013 1 0.711 -1.38 -1.38 0 0
#> # … with 127,730 more rows, 1 more variable: model <chr>, and abbreviated
#> # variable names ¹​quarter_f, ²​depth_sc, ³​est_non_rf
pred_num2 <- pred_num #%>%
#filter(year %in% c(2014, 2016, 2018, 2020, 2022))
plot_map_west +
geom_raster(data = pred_num2, aes(X, Y, fill = exp(est))) +
scale_fill_viridis_c(
trans = "sqrt",
name = "Count",
# trim extreme high values to make spatial variation more visible
na.value = "yellow", limits = c(0, quantile(exp(pred_num2$est), 0.99))
) +
facet_grid(model~year) +
ggtitle("Prediction (fixed effects + all random effects)",
subtitle = paste("Maximum estimated biomass density =", round(max(exp(pred_num2$est)), digits = 3))) +
geom_sf(size = 0.1)
ggsave("figures/poisson_maps.pdf", units = "cm", width = 20, height = 20)
# Index
dd <- num %>%
group_by(model, year) %>%
summarise(mean_num = mean(number))
#> group_by: 2 grouping variables (model, year)
#> summarise: now 30 rows and 3 columns, one group variable remaining (model)
p1 <- sims_num %>%
filter(.iteration < 26) %>%
ggplot() +
facet_wrap(~model, ncol = 3, scales = "free") +
geom_ribbon(data = sim_num, aes(year, est, ymin = lwr, ymax = upr), alpha = 0.2) +
geom_line(aes(year, .value, group = .iteration), alpha = 0.2) +
geom_line(data = dd, aes(year, mean_num), color = pal[8], size = 0.75, linetype = 2) + # add data
geom_line(data = sim_num, aes(year, est), size = 0.75) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 5)) +
theme_plot() +
theme(plot.margin = unit(c(0, 0.4, 0, 0.4), "cm")) +
xlab('Year') +
ylab('Mean count per tow') +
NULL
#> filter: removed 14,250 rows (95%), 750 rows remaining
year_diff <- sims_num %>%
filter(year %in% c(min(year), max(year))) %>%
pivot_wider(names_from = year, values_from = .value) %>%
mutate("2022-2013" = `2022` - `2013`) %>%
pivot_longer(c(`2013`, `2022`, `2022-2013`)) %>%
ungroup()
#> filter: removed 12,000 rows (80%), 3,000 rows remaining
#> pivot_wider: reorganized (year, .value) into (2013, 2022) [was 3000x4, now 1500x4]
#> mutate: new variable '2022-2013' (double) with 1,500 unique values and 0% NA
#> pivot_longer: reorganized (2013, 2022, 2022-2013) into (name, value) [was 1500x5, now 4500x4]
#> ungroup: no grouping variables
diff_sum <- year_diff %>%
filter(name == "2022-2013") %>%
group_by(model) %>%
summarise(diff_lwr = round(quantile(value, probs = 0.025), digits = 2),
diff_med = round(quantile(value, probs = 0.5), digits = 2),
diff_upr = round(quantile(value, probs = 0.975), digits = 2)) %>%
mutate(label = paste("2022-2013 = ", diff_med, " [", diff_lwr, ",", diff_upr, "]", sep = ""))
#> filter: removed 3,000 rows (67%), 1,500 rows remaining
#> group_by: one grouping variable (model)
#> summarise: now 3 rows and 4 columns, ungrouped
#> mutate: new variable 'label' (character) with 3 unique values and 0% NA
# Left join so that we can color "sig"
diff_sum <- diff_sum %>% mutate(sig = ifelse(diff_lwr > 0 & diff_upr > 0, "sig.", "not sig."),
sig = ifelse(diff_lwr < 0 & diff_upr < 0, "sig.", sig))
#> mutate: new variable 'sig' (character) with 2 unique values and 0% NA
year_diff <- left_join(year_diff, diff_sum)
#> Joining, by = "model"
#> left_join: added 5 columns (diff_lwr, diff_med, diff_upr, label, sig)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 4,500
#> > =======
#> > rows total 4,500
p2 <- ggplot(data = filter(year_diff, !name == "2022-2013"), aes(name, value)) +
geom_violin(alpha = 0.5, fill = "grey30", color = NA) +
geom_jitter(height = 0, width = 0.05, alpha = 0.2, color = "grey30", size = 0.1) +
geom_boxplot(fill = NA, width = 0.1, size = 0.5, color = "grey5",
outlier.shape = NA) +
theme_plot() +
theme(plot.margin = unit(c(0, 0.4, 0, 0.4), "cm")) +
facet_wrap(~ model, scales = "free", ncol = 3) +
xlab('Year') +
ylab('Mean count per tow') +
NULL
#> filter: removed 1,500 rows (33%), 3,000 rows remaining
p3 <- ggplot(data = filter(year_diff, name == "2022-2013"), aes(name, value, fill = sig)) +
geom_violin(alpha = 0.4, color = NA) +
geom_jitter(height = 0, width = 0.05, alpha = 0.3, color = "grey30", size = 0.2) +
geom_boxplot(fill = NA, width = 0.1, size = 0.5, color = "grey5",
outlier.shape = NA) +
scale_fill_manual(values = pal[c(2, 6)], name = "") +
theme_plot() +
theme(plot.margin = unit(c(0, 0.4, 0, 0.4), "cm")) +
facet_wrap(~ model, scales = "free", ncol = 3) +
geom_text(data = diff_sum, aes(label = label), inherit.aes = FALSE,
x = -Inf, y = Inf, size = 2.5, hjust = -0.05, vjust = 2) +
xlab('') +
ylab('Diff. in count per tow') +
geom_hline(yintercept = 0, linetype = 2, size = 0.4) +
NULL
#> filter: removed 3,000 rows (67%), 1,500 rows remaining
p1 / p2 / p3 + plot_annotation(tag_levels = "A")
ggsave("figures/poisson_index.pdf", units = "cm", width = 20, height = 20)Biomass predictions
# Space
pred_biomass
#> # A tibble: 127,740 × 12
#> X Y year depth year_f quart…¹ depth…² est est_n…³ est_rf omega_s
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 3.15e5 6.25e6 2013 228 2013 1 0.693 -2.60 -2.60 0 0
#> 2 3.17e5 6.25e6 2013 228 2013 1 0.693 -2.60 -2.60 0 0
#> 3 3.19e5 6.25e6 2013 227 2013 1 0.675 -2.60 -2.60 0 0
#> 4 3.21e5 6.25e6 2013 228. 2013 1 0.702 -2.59 -2.59 0 0
#> 5 3.23e5 6.25e6 2013 226 2013 1 0.656 -2.60 -2.60 0 0
#> 6 3.25e5 6.25e6 2013 229. 2013 1 0.717 -2.59 -2.59 0 0
#> 7 3.27e5 6.25e6 2013 229. 2013 1 0.711 -2.59 -2.59 0 0
#> 8 3.29e5 6.25e6 2013 228. 2013 1 0.686 -2.60 -2.60 0 0
#> 9 3.31e5 6.25e6 2013 229 2013 1 0.711 -2.59 -2.59 0 0
#> 10 3.33e5 6.25e6 2013 229 2013 1 0.711 -2.59 -2.59 0 0
#> # … with 127,730 more rows, 1 more variable: model <chr>, and abbreviated
#> # variable names ¹​quarter_f, ²​depth_sc, ³​est_non_rf
pred_biomass2 <- pred_biomass #%>%
#filter(year %in% c(2014, 2016, 2018, 2020, 2022))
plot_map_west +
geom_raster(data = pred_biomass, aes(X, Y, fill = exp(est))) +
scale_fill_viridis_c(
trans = "sqrt",
name = "Biomass density [kg]",
# trim extreme high values to make spatial variation more visible
na.value = "yellow", limits = c(0, quantile(exp(pred_biomass2$est), 0.99))
) +
facet_grid(model~year) +
ggtitle("Prediction (fixed effects + all random effects)",
subtitle = paste("Maximum estimated biomass density =", round(max(exp(pred_biomass$est))))) +
geom_sf(size = 0.1)
ggsave("figures/tweedie_maps.pdf", units = "cm", width = 20, height = 20)
# Index
p1 <- sims_biomass %>%
filter(.iteration < 26) %>%
ggplot() +
facet_wrap(~model, ncol = 3, scales = "free") +
geom_ribbon(data = sim_biomass, aes(year, est, ymin = lwr, ymax = upr), alpha = 0.2) +
geom_line(aes(year, .value, group = .iteration), alpha = 0.2) +
geom_line(data = sim_biomass, aes(year, est), size = 0.75) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 5)) +
theme_plot() +
theme(plot.margin = unit(c(0, 0.4, 0, 0.4), "cm")) +
xlab('Year') +
ylab('Biomass [kg]') +
NULL
#> filter: removed 14,250 rows (95%), 750 rows remaining
p1
year_diff <- sims_biomass %>%
filter(year %in% c(min(year), max(year))) %>%
pivot_wider(names_from = year, values_from = .value) %>%
mutate("2022-2013" = `2022` - `2013`) %>%
pivot_longer(c(`2013`, `2022`, `2022-2013`)) %>%
ungroup()
#> filter: removed 12,000 rows (80%), 3,000 rows remaining
#> pivot_wider: reorganized (year, .value) into (2013, 2022) [was 3000x4, now 1500x4]
#> mutate: new variable '2022-2013' (double) with 1,500 unique values and 0% NA
#> pivot_longer: reorganized (2013, 2022, 2022-2013) into (name, value) [was 1500x5, now 4500x4]
#> ungroup: no grouping variables
diff_sum <- year_diff %>%
filter(name == "2022-2013") %>%
group_by(model) %>%
summarise(diff_lwr = round(quantile(value, probs = 0.025), digits = 2),
diff_med = round(quantile(value, probs = 0.5), digits = 2),
diff_upr = round(quantile(value, probs = 0.975), digits = 2)) %>%
mutate(label = paste("2022-2013 = ", diff_med, " [", diff_lwr, ",", diff_upr, "]", sep = ""))
#> filter: removed 3,000 rows (67%), 1,500 rows remaining
#> group_by: one grouping variable (model)
#> summarise: now 3 rows and 4 columns, ungrouped
#> mutate: new variable 'label' (character) with 3 unique values and 0% NA
# Left join so that we can color "sig"
diff_sum <- diff_sum %>% mutate(sig = ifelse(diff_lwr > 0 & diff_upr > 0, "sig.", "not sig."),
sig = ifelse(diff_lwr < 0 & diff_upr < 0, "sig.", sig))
#> mutate: new variable 'sig' (character) with one unique value and 0% NA
year_diff <- left_join(year_diff, diff_sum)
#> Joining, by = "model"
#> left_join: added 5 columns (diff_lwr, diff_med, diff_upr, label, sig)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 4,500
#> > =======
#> > rows total 4,500
p2 <- ggplot(data = filter(year_diff, !name == "2022-2013"), aes(name, value)) +
geom_violin(alpha = 0.5, fill = "grey30", color = NA) +
geom_jitter(height = 0, width = 0.05, alpha = 0.2, color = "grey30", size = 0.1) +
geom_boxplot(fill = NA, width = 0.1, size = 0.5, color = "grey5",
outlier.shape = NA) +
theme_plot() +
theme(plot.margin = unit(c(0, 0.4, 0, 0.4), "cm")) +
facet_wrap(~ model, scales = "free", ncol = 3) +
xlab('Year') +
ylab('Biomass [kg]') +
NULL
#> filter: removed 1,500 rows (33%), 3,000 rows remaining
p3 <- ggplot(data = filter(year_diff, name == "2022-2013"), aes(name, value, fill = sig)) +
geom_violin(alpha = 0.4, color = NA) +
geom_jitter(height = 0, width = 0.05, alpha = 0.3, color = "grey30", size = 0.2) +
geom_boxplot(fill = NA, width = 0.1, size = 0.5, color = "grey5",
outlier.shape = NA) +
scale_fill_manual(values = pal[c(6, 2)], name = "") +
theme_plot() +
theme(plot.margin = unit(c(0, 0.4, 0, 0.4), "cm")) +
facet_wrap(~ model, scales = "free", ncol = 3) +
geom_text(data = diff_sum, aes(label = label), inherit.aes = FALSE,
x = -Inf, y = Inf, size = 2.5, hjust = -0.05, vjust = 2) +
xlab('') +
ylab('Diff. biomass') +
geom_hline(yintercept = 0, linetype = 2, size = 0.4) +
NULL
#> filter: removed 3,000 rows (67%), 1,500 rows remaining
p1 / p2 / p3 + plot_annotation(tag_levels = "A")
ggsave("figures/tweedie_index.pdf", units = "cm", width = 20, height = 20)knitr::knit_exit
#> function (append, fully = TRUE)
#> {
#> if (missing(append))
#> append = if (out_format(c("latex", "sweave", "listings")))
#> "\\end{document}"
#> else if (out_format("html"))
#> "</body>\n</html>"
#> else ""
#> .knitEnv$terminate = append
#> .knitEnv$terminate_fully = fully
#> invisible()
#> }
#> <bytecode: 0x7faec3c332a8>
#> <environment: namespace:knitr>